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GALERKIN FINITE ELEMENT METHOD FOR GENERALIZED 
FORCHHEIMER EQUATION OF SLIGHTLY COMPRESSIBLE FLUIDS IN 

POROUS MEDIA 

THINH KIEU t 

Abstract. We consider the generalized Forchheimer flows for slightly compressible fluids. Using Muskat’s 
and Ward’s general form of Forchheimer equations, we describe the fluid dynamics by a nonlinear degenerate 
parabolic equation for the density. We study Galerkin finite elements method for the initial boundary value 
problem. The existence and uniqueness of the approximation are proved. The prior estimates for the solutions 
in L°°(0, T; L q (fl)), q > 2, time derivative in L°°(0, T; L 2 (£l)) and gradient in L°°(0, T; VF 1,2_a (Q)), with 
a E (0,1) are established. Error estimates for the density variable are derived in several norms for both 
continuous and discrete time procedures. Numerical experiments using backward Euler scheme confirm the 
theoretical analysis regarding convergence rates. 
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1. Introduction . Fluid flow through porous media has been studied in many fields such 
as chemistry, physics, engineering, and geology. The most common equation to describe fluid 
flows in porous media is the Darcy law 


-Vp=-v, (1.1) 

K 

where p, v, p, k are, respectively (resp.), the pressure, velocity, absolute viscosity and permeability. 

For the Reynolds number large, for instance the high-velocity flow ©ESS, Darcy’s law 
is not valid any more. A nonlinear relationship between velocity and gradient of pressure 
is introduced, as suggested by Forchheimer in [TblUfi] . by adding the higher order term of 
velocity to the Darcy law. It is known as generalized Forchheimer laws, see [5] l30|I5Ill34|l36j and 
references therein. Forchheimer established the following three nonlinear empirical models: 

— Vp = av + b\v\v, —X7p = av + b\v\v + c\v\ 2 v, — Vp = av + d\v\ m ~ 1 v, m € (1,2). (1.2) 

Above, the positive constants a, 6, c, d are obtained from experiments. 

The generalized Forchheimer equation of ( 0 ) and m were proposed in mnm of the 
form 

N 

- Vp = y^ j q i \v\ ai v. (1.3) 

i=0 

These equations are analyzed numerically in Bilging, theoretically in mmrnm for 
single phase flows, and also in [221125] for two phase flows. 

In order to take into account the the dependence on density in generalized Forchheimer 
equation, we modify m using dimension analysis by Muskat J5U] and Ward [35]. They 
proposed the following equation for both laminar and turbulent flows in porous media: 

— Vp = F(v a K—p a ~ x where F is a function of one variable. (1.4) 

In particular, when a = 1,2, Ward (3S3 established from experimental data that 

— Vp = — v + cp— y=\v\v, where cf > 0. (1.5) 
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Combining m with the suggestive form (hd for the dependence on p and v, we propose 
the following equation 


N 

-Vp = '^2,a i p oli \v\ ai v, (1.6) 

i=0 

where N > 1, op = 0 < a.\ < ... < Un are real numbers, the coefficients ao, are 

positive. Here, the viscosity and permeability are considered constant and we do not specify 
the dependence of a^s on them. 

Multiplying both sides of previous equation to p, we obtain 

g(\pv\)pv = -pVp, (1.7) 


where the function g is a generalized polynomial with non-negative coefficients. More precisely, 
the function g : M + —> R + is of the form 

g(s) = aos“° + ais ai H- + aNS aN , s > 0, (1.8) 


where iV>l,ao = 0<a;i<---<Q ! iv are fixed real numbers, the coefficients ao,, a at are 
non-negative numbers with ao > 0 and a^v > 0. 

The equation of state which, for slightly compressible fluids, is 

l dp 1 

—— = — = const. > 0. 

p dp k 


Hence 


Vp=— pVp, or p\7p = tiS7p. 

Combining (11.81) and (11.91) implies that 

g(\pv\)pv = -nVp. 


(1.9) 


( 1 . 10 ) 


By rescaling coefficients of g we have 


g{\pv\)pv = -Vp. 


Hence 


where the function K : R + —> 


pv = -A'(|Vp|)Vp, 
is defined for £ > 0 by 


( 1 . 11 ) 


K(0 = 


1 


g(s{0) 

The continuity equation is 


, with s = s(£) being the unique non-negative solution of sg(s ) = £. (1-12) 


(j)pt + div(pv) = f. (1-13) 

where constant cf> £ (0,1) is the porosity, / is external mass flow rate . 

By combining (11.111) and (11.131) we obtain 

(j>Pt — V • (A'(|Vp|)Vp) = /. (1.14) 

Then by scaling the time variable in (11.141) implies 


Pt - V ■ (A(|Vp|)Vp) = /. 


(1.15) 
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The system of equations describing the fluid motion reduces to a scalar equation of density 
function. This is a nonlinear parabolic equation which degenerates as the gradient of density 
approach to infinity. 

For the existence and regularity theory of degenerate parabolic equations, see e.g. |7lH9l[27] . 
The popular numerical method for modeling flow in porous media are mixed finite element 
method. Arbogast, Wheeler and Zhang PQ first analyzed mixed finite element approximations 
of degenerate parabolic equation arising in flow in porous media. Not so long later, Woodward 
and Dawson in )37| study of expanded mixed finite element methods for a nonlinear parabolic 
equation modeling flow into variably saturated porous media. In recently years, the Galerkin 
finite element method for a coupled nonlinear degenerate system of advection-diffusion equations 
has studied in [MI- In their analysis, the Kirchhoff transformation is used to move the 
nonlinearity from coefficient K to the gradient and thus simplifies analysis of the equations. 
This transformation does not applicable for the equation (11.151) . 

In this paper, we combine techniques in EZH2DEH] utilizing the special structures of 
equation to obtain the error estimates for the solutions in several norms of interest. We obtain 
these results without any extra regularity assumption, though the order of error estimates are 
far from optimal order. 

The paper is organized as follows: In fJ5]we introduce notations, relevant results in mnm 
and suitable trace estimates in Lemma 12.41 for the nature of our equation. In ij3] we consider 
the semidiscrete finite element Galerkin approximation and an implicit backward difference 
time discretization to solve problem m ■ In we establish many bounds for solutions, its 
time derivative and gradient to problem (13.21) and (13.41) in Lebesgue norms. In ij5]we analyze 
two version of a Galerkin finite element approximations, the continuous Galerkin method and 
the discrete Galerkin method. Using the monotonicity properties of Forchheimer equation and 
the boundedness of solutions, the priori error estimates are derived for solution in L 9 -norm 
with 2 < q < oo and for gradient of solution in L 2_a -norm. Finally, in ij6] numerical examples 
using the Lagrange elements of order 2 are carried out in two-dimension. The results strongly 
support our theoretical analysis regarding convergence rates. 

2. Notations and preliminary results. Suppose that Cl is an open, bounded subset of 
R d , with d = 2,3,..., and has ^-boundary T = dCl. Let L 2 (Cl) be the set of square integrable 
functions on and (L 2 (Cl)) d the space of d-dimensional vectors which have all components in 
L 2 (Ct). We denote (•,•) the inner product in either L 2 ( f2) or (L 2 (Cl)) d that is 



The notation (•,•) will be used for the L 2 (dCl) inner-product and ||u|| iP = ||'it|| iP (Q) for 
standard Lebesgue norm of the measurable function. The notation ||-|| will means scalar norm 
IHIl 2 (Q) or vector norm |H| (L2(n))d . We denote \\u\\ LP{Lq) = |M| LP(0 ,T;L«(n)), 1 < P,q < oo 
means the mixed Lebesgue norm for a function u while = ||it|| iP , 0 , 1 — 

p,q < oo stand for the mixed Sobolev-Lebesgue norm of a function u. 

For 1 < q < +oo and m any nonnegative integer, let 


W m ’ q (Cl) = {it G L q (Cl),D q u G L q (Cl), | 9 | < to} 


denote a Sobolev space endowed with the norm 



| q;| <m 


Define U m (U) = W m ’ 2 (D) with the norm ||-|| m = ||-|| m 2 . 

Throughout this paper, we use short hand notations, 

A = (0, r), \\p(t)\\ = \\p(-,t)\\ L2{n) yt€l and p°(-) = p(-,0). 
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Also our calculations frequently use the following exponents 


and 


a-N _ deg(g) 
ctN + 1 deg(g) + 1 ’ 


0 = 2 — a, 


2 - 




1 - 


0 a a 

0-1’ 7= 2-a = 0' 


( 2 . 1 ) 


( 2 . 2 ) 


The arguments C, Co, C*i,... will represent for positive generic constants and their values 
depend on exponents, coefficients of polynomial g, the spatial dimension d and domain Q, 
independent of the initial and boundary data, size of mesh and time step. These constants 
may be different place by place. 

Lemma 2.1 (cf. [211171 ). The function AT(0 has the following properties 

(i) K : [0,oo) —> (0, a-o” 1 ] and it decreases in £, 

(ii) For any n > 1, the function A"(OO 1 increasing and AT(0£" > 0 

(Hi) Type of degeneracy 


(i + C)° ~ ~ (i + 0“’ 

(2.3) 

(iv) For all n > 1, 


c 3 (c~ a - 1) < k( or < car* -0 , 

(2.4) 

(iv) Relation with its derivative 


- aK( 0 < K\ Of < 0, 

(2.5) 


where ci, 02,03 are positive constants depending on O. and g. 

We define 

£ 2 

if(£) = f K(\fs)dx, for £ > 0. (2.6) 

Jo 

The function iJ(0 can compare with £ and K(tf) by 

Km 2 < mo < mm 2 - ( 2 . 7 ) 

For the monotonicity and continuity of the differential operator in (11.1511 we have the 
following results. 

Lemma 2.2 (cf. |17l). One has 

(i) For all y, y' G M d , 

( K (\y'\)y' - K (\y\)y) ■ (v' -y)>(P- l)if(max{|g|, \y'\})\y' - y \ 2 ■ (2.8) 

(ii) For the vector functions Si,S 2 , there is a positive constant C such that 

(AT(|si|)si - A'(|s 2 |)s 2 ,si - s 2 ) > Cw||si - s 2 ||^ ( n) , ( 2 -9) 

where u = (l + max{||s 1 || £ ^ (n) , ||s 2 ||i/3(n)}) “• 

Lemma 2.3 (cf. [26] 1 . For all y,y' G R d . There exist a positive constant C depending on 
polynomial g, the spatial dimension d and domain S2 such that 

\K(\y'\)y' -K(\y\)y\<C\y' ~y\- (2.10) 

Next we derive trace estimates suitable to our nonlinear problem. 
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Lemma 2.4. (i) Let q > 2. Assume v(x) is a function defined on Cl. If \v\ q 1 £ W 1 ’ 1 ^) 
then for all e > 0, 

Jr M q " X dxr < C IMI^n) + 5 J Q \ v \ q ~ 2 \ Vv \ 0dx + Ce~i^ IMI^n) • ( 2 - n ) 


(ii) IfuE L°°(r) and |v| G W 1,1 (n) then for all e > 0, 


'I Lf>(Q) 


c 


( £ 1 IMlL»(r)+ e p 1 ll u llL~(r) 


|{u,w}| < £ (\\v\\ 2 + ||Vv|| 

In particular case, 

IMI< J(lbll 2 + liv^||^ (n) )+c(i + || u || 

Above, C is a positive constant independent of u, v, e. 

Proof. We recall the trace theorem 


)■ 


l°°( r) ■ 


( 2 . 12 ) 


(2.13) 


\4>\dx < C [ \(f\dx + C [ \V<j)\dx, 

J n Jq 


for all 4> £ W 1,1 (f2), where C are positive constants depending on Cl. Applying the Trace 
theorem to 4> = H ?_1 shows that 


[ \v\ q ~ 1 da < C [ \v\ q ~ x dx + C [ \v\ q ~ 2 \S7v\dx. 
Jr J n J n 


(2.14) 


(2.15) 


Note that ^ + j = 1 and j = Using Young’s inequality with exponent /3 and A, we find 
that for £ > 0 

\v\ q ~ 2 \C7v\ =£?\p h \ !L r\\'v\-E-Ti\v\ 3 ^ 

< £\v\ q - 2 \S7vf +C£-^\v\ q ~ 2 . 

Combining (12. 1411 . (12.151) we obtain 

f \v\ q ~ 1 da<C f \v\ q ~ 1 dx + e f \v\ q ~ 2 \Vv\^dx + Ce~'^ rT f \v\ q ~ 2 dx. (2-16) 

Jr i/ >/ fj >/ 

Inequality (12.111) follows by using Holder’s inequality to the first term of the right hand 
side in (12.161) . 

(ii) We have 


|(zt,w)| < |MU (r) J \v\da. 

For all <5 > 0, using (12.161) with q = 2 gives, 

f \v\da < C f \v\dx + 5 || Vu||^ 3 + CS~ p- 1 . 
Jr J a 

Applying Young’s inequality leads to 


(2.17) 


\v\da < 5 |M|" + C8~ L + 5 ||Vv||^ + C5~fc, 


which gives 


\M\ < ll«llL- ( r) («J|H| 2 + C'<r 1 + <J||Vt;||^ +C8-T&). 


If \\u\\ LOO[r) = 0 then (12.121) clearly holds true. 

Otherwise, selecting 5 = e IMI^co/p) and the fact that + 1 = A we obtain estimate 

(EH. 

Estimate (12.131) follows by choosing £ = 1/4 in (12.121) and using Young’s inequality. □ 
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3. The Galerkin finite element method. Our aim is to study equation (11.151) for 
density of slightly compressible fluids in bounded domain in porous media. The fluid flows are 
subject to some conditions on the boundary. 

We consider the initial boundary value problem associated with (11.151) , 

(p t - V • (K(\Vp\)S7p) = f in Ox/, 

< p(x,0) = p°(x) in O, (3.1) 

[it (\X7 p\)X7 p • v + ip = 0 on r x/, 

where p°(x) and ili(x,t) are given initial and boundary data, respectively. 

The variational formulation of m is defined as the following: 

Find p : [0, T\ —> W = i/ 1 (0) such that 

{pt,w) + (iF(|Vp|)Vp, Vw) = -(ip, w) + w £ W (3.2) 

with p(x, 0) = p°(x). 

Let {Th}h be a family of quasi-uniform triangulations of O with h being the maximum 
diameter of the element. Let Wh the space of discontinuous piecewise polynomials of degree 
r > 0 over Th. It is frequently valuable to decompose the analysis of the convergence of finite 
element methods by passing through a projection of the solution of the differential problem 
into the finite element space. Here we use the standard L 2 -projection operator (see [6]) 7r : 
W —> Wh, satisfying 


(-KW,v h ) = (w,v h ), \/weW,v h &W h - 

These projections have well-known approximation properties (see [5J[25]). 

(i) ||7ru;|| < ||w|| holds for all w € L 2 (fi). 

(ii) There exist positive constants C \, 6b such that 

ll™ ^ w\\l<(Q) < Cih™ |Mlvr«>,(3-3) 

for all w £ W m,q (Q), 0<m<r + l,l<q<oo. Notation ||-|| m denotes a standard norm in 
Sobolev space W m ’ q (fl). In short hand, when q = 2 we write (13.31) as 

\\irw - w\\ < Cih,™ \\w\\ m and \\nw - w|| Loo(n) < C 2 h m |ML +1 ■ 

Replacing the original density by its approximation, the semidiscrete formulation of (13.21) 
can read as following: Find pn ■ [0,T] —» Wh such that 

(ph,t, w h ) + (K(\\7p h \)\7ph, ^w h ) = ~{i>, w h ) + (/, w h ), w h £ W h (3.4) 

with initial data = 7rp°(x). 

Let N be the positive integer, t 0 = 0 < ti < ... < fjv = T be partition interval [0,T] 
of N sub-intervals, and let At = t n — t n -\ = T/N be the n-th time step size, t n = nAt and 
ip n = ip(-,t n ). The discrete time Galerkin finite element approximation to (13.21) is defined as 
follows: 

Find p£ Wh, n = 1, 2,..., N, such that 

( P * A f ,v>h) + {K(\Vpl\)Vpl,Vw h ) = ~(r,w h ) + (f n ,w h ), Vw h £ W h , (3.5) 


with initial data are chosen as follows: p° h (x) = np°(x). 
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4. A priori estimate for solutions. We study the (13.21) . (13.41) and (13.51) equations for 
the density with fixed the functions g(s) in (11.71) and (11.81) . Therefore, the exponents a, and 
coefficients a; are all fixed, and so are the function K (£), 7J(£) in (11.121) . (12.61) . 

With properties m , ED , ED, the nronotonicity (12.81) . and by classical theory of 
monotone onerators (28U331I38I. the authors in m proved the global existence and uniqueness 
of weak solution of equation (13.21) . Moreover p £ <7(0, T, L q (Vl)) D if oc (0, T, W 1 ’?^)), q > 1 
and p t £ L^ c (0, T, (W 1,/3 (f2))') D Lf oc (0, T, L 2 (fl)) provided the initial, boundary data and 
/ sufficiently smooth. For a priori estimate, we assume throughout this paper that 
belong to C([0, T], L°°(r)), p° £ L°°(fl) and / £ C([0, T], L°°(fl)) and the weak solution are 
sufficiently regularities both in x and t variables such that our calculations can be performed 
legitimately. 

The local existence of of approximate solution of (13.41) follows from Peano’s theorem. The 
global existence of ph relies on the known theory of ordinary differential equation and stability 
estimate which shows as follows. 

Theorem 4.1. Let q > 2 and ph be a solution of ED- Then there exist positive constant 
C such that 


\\Pkh 






c 


jfo 


qA 

2 

L°°(r) 



(4.1) 


Proof. Selecting Wh = \ph\ q 1 sign(p ft ) in (13.41) . we obtain 

\j t WphWI, = ~{q - 1) (A-(|Vp,|)|Vp,| 2 , \ Ph \ q ~ 2 ) 

- (V 1 , 1 Pa. 1 1 sign (p/i)) + (/, |p^| 9_1 sign(p ft )) . 

It follows from (12.41) and Holder’s inequality that 

-(9-1) (iqiV^DIV^I 2 , \p h \ q ~ 2 ) < -Co [ {\S7ph\^ — l)\ph\ q ~ 2 dx 

Jo. 

= -c 0 [ \^Jph\j 3 \ph\ q ~ 2 dx + c 0 f \p h \ q ~ 2 dx (4.3) 
Jn J n 

<-c 0 [ \Vp h f\p h \ q - 2 dx + C\\p h \\ q - 2 , 

J n 

where co = Cs(q — 1). 

According to (12.111) in Lemma 12.41 


\-(^,\Ph\ q 1 sign(p /l ))| < ||V>|| I ,=o (r) f \ph\ q 1 da 

j r 

< { \\Ph\\ q L^ + e / \Ph\ q ~ 2 \^Phfdx + £~-^ \\p h \\ q ~ 2 y (4.4) 

J 

Using Young’s inequalities with exponent q/(q — 1) and q , we obtain 

(/, |p/ l | 9_ ‘ 1 sign(ph)) < C \\phWl- 1 + C \\f\\ q Lq . (4.5) 

Combining (14.31) . (14.41) . (14.51) and (14.21) gives 


-^l\\Ph\\ q Lq <-ca I \\7 Ph f\ Ph \ q 2 dx + C\\p h \\ q Lq 2 
q at Jn 

+ C\\fj\\ L o. ir) [\\Ph\\ q L ~ q 1 +£ f \p h \ q - 2 \S7p h \ p dx + e-T^ \\p h \\ q ~ 2 } 

J £1 

+ C\\p h \\ q L - 1 + C ||/|| 9 (4.6) 
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The case || V'(^) II(r) = by the means of Young’s inequality in (14.61) we obtain 

^ \\Ph\\ q Lq + qco WPhf\Ph\ q ~ 2 dx < C \\p h \\ q Lq + C{1 + \\f\\ q Lq ). 
Applying Gronwall’s inequality to above differential inequality implies that. 


(4.7) 


\Ph\\L°°(I,Li(Q)) + d c 0 


0 C(t-T) 


o J n 


\^Ph\ P \ Ph \ q - 2 dxdT 


< c U\\Uu)+ c 


jf(‘ + 


l|«(fi) 


(4.8) 


dt. 


Consider HV’C^)lli°°(D ^ 0- Selecting £ = ll'^lli^qr) i n CE2b we obtain 

\\ph\\h<~ [ \s/p h \ p \p h \ q - 2 dx+c\\p h r L - 2 

z Jn 

{ IIp^III, + IIV’lli?(r) Uphill® } + C IIp^II^o + c | 


d_ 

dt 


c 


n 

l °°( r) 


(4.9) 


I Li ■ 


By Young’s inequality, 

j t \\ph\\ q Lq J Q \Vphf\Ph\ q ~ 2 dx 

<C\\p h \\l q +c(l + \mi~<r) 
<c|K||l, + c(i + M^ 

Solving this differential inequality shows that, 


Q0 

I 20-1) 

lL°°(r) 


II. 


(4.10) 


1 1 


)• 


\ph\\ q L , 


+ H0 e C(t-r)\S7p h \P\p h \1-* dxdT 

z Jo Jn 


r l / ql3 

— Ik° III, + C J (l + Ill’ll L°°(r) + 

It is easy to see that in both above cases 

||P/l|lI? < \\Ph W^q + C j (l + 


O 

20 - 1 ) 

L°°(r) 


III, 


III, 


y 'jdr . 


^jdr. 


(4.11) 


Using inequality (a + &) 1 / 9 < 2 1 / q (a 1 ^ q + b 1 / q ) with 


II on? 

a = IKIIl, 


and b = 


/ ( 1 + W 


O 

20 - 1 ) 

i~(r) 


II, 


dr. 


we obtain 

IIP/t|lt°°(/,L,(n)) 

Note that 


SCIKIUnj+C 


jf(‘ + 


qA 

li~(r) 


l|,(n) 


dt 


\ph\\ = |kp°|| < ||p°||- 


(4.12) 


(4.13) 


Thus inequality (HU) holds. We finish the proof. □ 

Next we give estimates in L 2 -norm following directly from (14.111) and (14.131) with q = 2. 
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that 


Lemma 4.2. Under assumption of Theorem There exist a positive constant C such 


\\Ph\\L°°(I,L 2 (n)) + II Vp fc || LP(i,LP(a.)) 


< C\\p° +C 


(l + ll^)ll^ (r) + ||/W|| 2 )dt 


(4.14) 


Since the equation (13.41) can interpret as the finite system of ordinary differential equations 
in the coefficients of ph with respect to basis of Wh- The stability estimate (14.121) suffices to 
establish existence of Ph{t) for all t G I. 

For the uniqueness of approximation solution, assume that for i = 1,2, pj^ £ Wh is the 
solution of (13.41) . Let ph = pj^ ~ P^ then 

(Ph,t,w h ) + {KQVpPnVpW - K(\Vp™\)Vp™, Vw h ) = 0. 

Choose Wh = Ph, implies 

\j t \\Ph\\ 2 + (^dVp^DVpW - A'dVpfDVpf.VpW - Vp' 2) ) =0. 

According to (12.91) we find that 


(iqiVp^DVpW 


K{\S7p™\)Vpl 


.Vp^ 



> c 


l + maxlHVpWll^JIVpflU,} ||Vp h || 


L« 


Hence 


\A 11/9/11,2 + C I 1 + maxlHVp^lU,, ||Vpi 2) || L3 }] “ \\VphWl? < 0. 


Integrating in time from 0 to t we have 


I \Phf+C / [(l + maxlllVp^lU^HVpf \\ Lf> }] ° || Vp h f L , dt < |M0)|| 2 = 0. 


/o L 


Thus 


Ph = 0, Vph = 0 V(x, t) € ft x I. 


Next, we find estimates for Vph- We define 

g(t) = 1 + ll'0(*)llioo ( r)» Ht) = 1 + \\Mt)\\L°°(r)i fc(t) = p(t) + ||/(t)|| 2 . (4.15) 

Theorem 4.3. Let Ph be a solution to the problem (Ell). Then there exist constant positive 
constant C satisfying 

\\Ph,t\\L 2 (I,L 2 (Q.)) + ll^frlli£(/,L£(n)) T \\^Ph\\L°°(I,L?<(n)) + ||P/i|li=o( / i 2(Q)) < CM. (4.16) 


where 


M = max gft) + [ [h(t ) + k(t)\dt + [ [ kMdrdt 

*e[o,T] J g J 0 J 0 

01 


+ Imp 


oil / 3 


I Lf>(Q) 


W 


lL 2 (r) 


li 2 (r) 


(4.17) 
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Proof. Choosing Wh = Ph,t in (13.41) leads to 
\ d f 

IIpmII 2 + 2 J H(x,t)dx = - ( ip,Ph,t } + (/,Pm) 

= WbPft) + + (f,Ph,t)- 

where H(x,t) = H(\\7ph{x,t )|) is defined in (12.61) . 

With q = 2 then from (14.101) . 

j t \\Ph \\ 2 + co || VphC < C \\p h \\ 2 + c(l + ||V-||^ (r) + ll/ll 2 ). 

Let 

£(t) = [ H(x,t)dx+ \\p h \\ 2 + 2 {if,p h ). 

Jn 

Adding two inequalities (14.181) and (14.191) gives 

\\Ph,t\\ + co IIVp^ll^ + — —£(t) < (ip t ,ph) + (/, Ph,t) + C Up*.|| +Ck(t). 
Using (12.131) and Cauchy’s inequality then 


(4.18) 


(4.19) 


\\ph,t\\ + co ||Vp fc ||^ + ~—£(t) < - ^||p?i|| + \\Vp h f Ll> ) + C ^1 + ||V ; i|lL“(r) 


+ O ll/ll + o \\ph,t\\ + C||P?i|| +Ck(t). 


Absorbing ||/|| 2 to k(t), integrating in time, we obtain 


[ \\Ph,t\\ 2 dt + c 0 [ \\VPh\f L fi dt + f H(x,t.)dx+\\p h \\ 2 <-2{i/;,p h ) 
J o Jo Jn 


+ C f || p h fdt + C f [h{t) + k(t)]dt + £{ 0). 
Jo Jo 


(4.20) 


Applying (12.121) to the first term of the left hand side in (14.201) and using the fact in (12.71) that 
co{\^7Ph\ d - 1) < H(x,t) < 2c 2 \\/p h \ 0 , we have 


IIpmII dt + c 0 HV^H^dt + csIlVpfcll^ + IKir 


L-(r) 


+ e f>- 


\ 

L°°(r) 


<2e{\\p h \\ 2 + \\V Ph f Lfi )+c{e~ 1 \m 

r T r T 

+ C || Ph \\ 2 dt + C [h(t) + k(t)]dt + £( 0). 

J 0 Jo 

Thus by taking e = min{c 3 , l}/4, and using Young’s inequality, 

\\Ph,tf dt + c 0 \\V Ph\\L?dt + ^\\V p^e + ^\\p h \\ 2 <Cg(t) 

r T r T 

+ C II p h fdt + c [h(t) + k(t)]dt + £( 0). 

Jo Jo 


Note from (14.111) with q = 2 that 


\\p h fdt<C[T\\p h (0)\\ 2 + 


ff 

JO Jo 
r T r t 


k(r)drdt 


< C T hr + 


k{r)drdt 


>o Jo 




(4.21) 


(4.22) 
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and from (12.71) . (12.41) that 

£(o) < c (||^f + ||v^||^ (n) + ||^(o)|| i2(r) K|| i2(r) ) 

< c (||p°|| 2 + ||Vp°||^ (n) + ||^(0)|| i2(r) ||p°|| i2(r) ). 

The left hand side of (14.221) is bounded by CM which implies (14.161) . □ 

Now we prove that the time derivative of density is also bounded. 

Theorem 4.4. Let 0 < to < T, ph be a solution to the semidiscrete problem (13.41) . Then 
there exist constant positive constant C such that 

\\ph,tfL~(i,v> m < Ct-'M+CM+C j\(t) (l + ||/ t (f)|| 2 ) dt. (4.23) 

where M and h(t ) are defined in (14.171) and (14.151) respectively. 

Proof. Differentiating d3~H) with respect t yields that 

(. Ph,tt,w h ) + {K(\Vp h \)Vph,t, Vw k ) 

= - ^'(1 Vpft|) Vp h , Vwh'j - (iftMh) + {ft,w h ). 

Choosing Wh = Ph,u we obtain 


Id 
2 dt 


\p h ,tf + K 1/2 {\S7p h \)S7p Ki 


= - ( ^(IVjOhl) V ^ V ^ M Vp fe , Vp fe , t ) + UuPh,t) - (A,Ph,t) ■ 


Using (12.5f) . 


Thus 


- ( K\\S7p h \) ^ ph ^ Kt ^p h yp h ^ 


< a 


A' 1/2 (|Vp,|)Vp,, t . (4.24) 


ll/MI 2 + (!-«) H^(|Vp,|)V P/l , t || 2 < ( ft,Ph,t ) - (i>t,Ph,t) ■ 
In virtue of Cauchy’s inequality, for all e > 0 

( ft,Ph,t ) <£\\ph,t\\ 2 + C£ ~ 1 Wftf ■ 

Using Trace Theorem we obtain, 

\(A,Ph,t)\ < ||^t|| L «. ( r) [(I Ph,t\A) + (|Vp M |,l)]. 

Again Cauchy’s inequality gives that for all e, £\ > 0 

(IpmM) < £ \\Ph,t \\ 2 +Ce~ 1 , 

and 

(|Vp M |, 1) < e 1 (K(\^7p h \)\S7p htt \ 2 , 1) + Cef\K~\\S7p h \, 1). 
By using (12.31) and (1 + x) a < 1 + x a , x > 0 imply 

(|Vp M |, 1) < £i (K(\Vp h \)\Vp h , t \ 2 , 1) + Cef\l + |Vp h |°, 1) 


< £i 


K 1/2 (\Vp h \)Vp h ,t + Cef\ 1 + ||V Ph ||“,). 


(4.25) 

(4.26) 

(4.27) 

(4.28) 


(4.29) 
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Combining (14.271) . (14.281) and (14.291) shows that 
\{i>t,Ph,t)\ < HV’tIIL°° ( r) { e \\Ph,t\\ 2 +Ce~ l 


+ £ 1 


K 1/2 (\Vp h \)Vph,tf + Cef(l + ||Vp h ||“,)}. 


(4.30) 


It follows from (14. 25f) . (14.261) and (14.301) that 
^||pmII 2 + (i-«)||^ 1 / 2 (|v^|)v^, " 2 


< mill/ 


yr) { £ IIpmII 2 + C £ 


+ £\ 


K 1/2 (\Vp h \)Vp h ,t\\ + Ce r'(l + IIVphH^)} + e ||p M f + Ce 


2 i r>^~ 1 ll rjl 2 


(4.31) 


If ll^t(*)ll L ~ ( r) = 0 then 

^\\ph,tf+ ( 1 ~ a ) K 1 / 2 (\Vp h \)Vp h ,t < eHpmII 2 + C , e " 1 

Dropping the second term of the left hand side in above estimate, selecting e = 1/2 and, 
for 0 < if < to < t < T integrating from t' to t give 


WphAM 2dT < IIpm(OII 2 + J f \\Ph,tf dT + G J ll/tll 2dT 

<IIpm(OH 2 + f \\Ph,t \\ 2 dT + C fwftfdr. 

Jo Jo 


(4.32) 


Now integrating in t' from 0 to to, 

to \\ph,t(t)\\ 2 dT < WPhA^f +t °{J o \\Ph,t\\ 2 dT + C J q Wftfdry (4.33) 

Using (14.161) . we find that 

*o WphAVf d,T<M+Ct 0 (M + J Wftfdry (4.34) 

which holds (14.231) . 

Consider ||-0tllz.°°(r) ^ 0. Selecting £l = ^ ||^||^i (r) and e = |(||i/’t|| L oo (r) + l)" 1 hr 
(14.311) yields 

j t \\PhA 2 + (i -<0 H^dVp^DVp^f < \\ph,t\\ 2 + c WMl-w \\vph\\h 

+ C (jlV’tllLoo( r ) + l) IIL°°( r ) + ll/*l| 2 ) + C ll^tl| 2 <»(r) 

< IIpmII 2 + W^PhWLd + cz{t), 

where 

Z{t) = (llV’i|lL«>(r) + l) (lIV 1 *IIL°°(r) + ll/*l| 2 ) + 11^*11 l°°( r) + HV^tllz)°°(r) ■ 

For t > to > t' > 0. Ignoring the the nonnegative term in the left hand side of above 
inequality, integrating from t' to t and then integrating in t' from 0 to to, we find that 

toWPhA 2 < [ \\Ph,ttf)\f dt '+to [ \\Ph,tW 2 + \\^Ph\f L p dr + Cto [ Z (r )dr. 


’ o L 


By virtue of (14.161) . 


/ o L 


IIpmU +IIV^II^ dt < CM, / \\ Pht t(t'W dt <CM. 
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Therefore 


to \\Ph,t\\ 2 <CM + CMt 0 + Ct 0 f Z{r)dr. 

Jo 

We estimate Z- term by 

Z(t) < (||V’t|li,°°(r) + l) ll/tll + C (l + ll^t IIL°o(r)) 

c(l + \\A\\U { r)) (l + ll/dl 2 )- 


< 


(4.35) 


(4.36) 


The inequality (14.2311 follows from (14.351) and (14.361) . The proof is complete. □ 

The following results can be proved by following the ideas of the proof of Theorem l4.ll 14.31 
and 14.41 

Theorem 4.5. Let 0 < to < T, and q > 2, p be a solution to problem m ■ Then there 
exist constant positive constant C such that 

pT x I 

IML-d.wCn)) <C\\p°\ \ LHn) +c( J Q [1 + mt-iD + \m\\Un)}dty, (4.37) 

WpA\l 2 (i,l 2 {q)) + \\^pf L s(i,Lf<(n)) + ll^ 7 PllL“(/,i/ 3 (n)) + l|p/iII(j,i 2 (r2)) - Oi'M, (4.38) 

pT 

\\pt\\ 2 L~> { i,L* m < Cto l M + CM + C h(t) (l + ||/ t (t)|| 2 ) dt, (4.39) 

J 0 

where h and A4 are defined as in (14.151) and (14.171) . 

5. Error estimates. In this section we will establish the error estimates between analytical 
solution and approximation solution in several norms. 

5.1. Error estimate for continuous Galerkin method. We will find the error bounds 
in the semidiscrete method by comparing the computed solution to the projections of the true 
solutions. To do this, we restrict the test functions in (13.21) to the finite dimensional spaces. 
Let 


X = P ~ Ph = (p ~ np) - ( p h - up) = d - O h . 


(5.1) 


Theorem 5.1. Let p,ph be two solution to problems m and dS3D respectively. Assume 
that p £ L°°(I, H r (fl) D W r ’P(Ll)), pt £ L 2 (/, H r (12)). Then there exist positive constants C 
independence of h such that 


\\p ~ P^IIl°o(/,L2(Q)) + M 2 || V(p — P/l)|| L 2(/ ii 3(Q)) < Ch' ||p|| L oo(/ iff r(f2)) 

+ CM 2 t> x h~ ||p||£l( W r,3(Q)) + CT 2 1T ||pt||i2(j )ff r(n)) , 

where Ai is defined as in (14.171) . 

Proof. Subtracting (13.41) from (13.21) we obtain the error equation 

{pt ~ Ph,uWh) + (^(|Vp|)Vp- K(\Vp h \)S7p h , S7w h ) = 0, w h £ W h . (5.3) 

Choosing Wh = Oh in (15.31) gives 

\j t \\0h\\ 2 + (K(\Vp\)Vp - K(\Vp h \)Vp h , Vp - V Ph ) 

= (d t ,0 h ) + (K(\Vp\)Vp - K{\Vp h \)Vp h , Vi?). 
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By the monotonicity of K(-) in (12.91) . there is a positive constant Co independence of h and t 
satisfying 

(if (|Vp|)Vp - K(\Vp h \)Vp h , Vp - \7 Ph ) > C 0 uj(t) ||Vx||^ - 
where u(t) = (! + || Vp||^ + || ^ph\\ L p)~ a ■ 

Observing from (14.161) that there is a positive constant C\ independence of h and t such 

that 


Thus 


w 1 (^) — (1 + ll^PlIiP + \\^ PhWLoY 1 

< ^ (i + \\vpf Lf >+ ii vphW^y < ci m\ 


(K (|Vp|)Vp - K(\Vp h \)Vp h , Vp - Wp h ) > WVxWle • 


(5.5) 


(5.6) 


By Cauchy’s inequality, for > 0 

Wt,O h )<Ce o'lltftf+eollM 2 , (5-7) 

Using (El, Holder inequality and then (14.161) we obtain 

(AT(|Vp|)Vp - if(|Vp h |)Vp fc , V0) < C(| Vp^- 1 + \Vp h f~\ |V0|) 

<C'(||Vp||^ + ||Vp h ||| <I )||Vt?|| i<I 

< C7(||Vpll^ + llVphll^)^ IIV^II^ 

< CM?* ||Vd|| i/3 . 

We obtain from (15.41) . (15.61) . (15.71) and (15.81) that 

~ ll^ll 2 + HVxIli, < CA 4 A ||W| L , + Ceq 1 ||tf t || 2 + £o \\ 0 h \\ 2 ■ 

Integrating in time from 0 to T and then taking sup-norm in time of resultant show that 


(5.8) 


IIM 


L°°(I,L*) 


+ M.-1 [ WVxWl'dt< CM** [ || W || L p dt + Ceq 1 f ||tf t || 2 dt + e 0 T \\6 h \\ 
Jo Jo Jo 


Selecting £o = we find that 


\\0 h \\U (I L2) + [ T \\VxWl, dt<CM& f T \\Vn LP dt + CT [ T \\M 2 dt. (5.9) 
Jo Jo Jo 

It follows directly from (I5.9|) that 

WOhh-^+M-i ||Vx|| ia(/ , ifl) < CA4 w ||W||| 1(/ ^ + CTi ||tf t || L2(J , i2) 

< CM^h^ ||p||J 1(JiWri(J) + CT^h r ||pt|| L2(W - 

(5.10) 


By triangle inequality, 


I L“(/,L 2 ) 


< M 


L°°(I,L*) 


m 


L°° (I,L 2 ) ■ 


This and (15.101) imply (15.21) . Thus the proof is complete. □ 
Now we give the error estimate in L q -norm for any q > 2. 


UT 2 ) • 
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Theorem 5.2. Let q £ (2,oo), p solve the problem (13.21) and ph solve the semidiscrete 
problem (13.41) . Assume that p £ L°°{I,W r ’ q {fl) D W r ’ 2q {Ll)), Vp £ L 2q {I, L 2q {fl)) , p t £ 
L q (I,W r,9 ( fi)). Then there exist constant positive constant C independence of h and T such 
that 


\\P Ph\\L°°(I,LQ(n)) — \\p 


7r PllL“(/,L'J(n)) + 5 II (p ' K P)t\\Li(I,Li(£l)) 

+ CA«T?-« ||V(p-7rp)|| i2s(J)i28(n)) , 


(5.11) 


where 


A = 


l + l|Vp||^(n) + IIV^||^ ( a) 


dt 


(5.12) 


Consequently, 


\\P Ph^L°°(I,Lt(Q)) — Ch \\p\\L^(I,W r -i(Q.)) + CT ih ||Pt|li< 1 (7 j VE r ^(n)) 


+ CA«T*—«h r 


(5.13) 


\L 2 i(I,W r ’ 2 i(Q)) ■ 


Proof. Choosing Wh = \9h\ q 1 sign0^ in (15.31) we have 

(t?t - e h , u l^l^^ign^) + {q- l)(K(\Vp\)Vp - K(\Vp h \)Vp h , \9 h \ q ~ 2 V9 h ) = 0. 

We rewrite as form 

(0m, l^l 9_1 sign^) + { q - i)(^(|Vp|)Vp- K(\Vp h \)Vp h , \e h \ q ~ 2 (y P - V Ph )) 

= {dt, Iflfcl^signflh) + {q- l){K{\Vp\)Vp - K{\S7 p h \)S7 p h , |6> h | 9 " 2 Vi?). 

The nronotonicity of K{-) in (12.81) provides 

(K{\Vp\)Vp - K(\Vp h \)Vp h , 1 9 h \*- 2 (Vp - V Ph )) >{P~ 1){K (OlVxl 2 , l^l 9 ” 2 ), (5.14) 

where £ = max{ | Vp|, \S7ph\}. Hence 


-4 II^IIL + (9 - m - 1)(A'(0|V X | 2 , l^l 9 " 2 ) 

qat 


< {dt, |flhl^signflfc) + {q- l)(K{\Vp\)Vp - K{\Vp h \)Vp h , \9 h \ q ~ 2 Vd). 


Using Young’s inequality with exponents q and 


{d t ,\9h\ q 1 sign 9 h )<q 


-l 


9-1 


1-9 


e l - q \\d t r Li +e\\9 h r Lq . 


Using (12.10|) . Holder’s inequality and Young’s inequality we obtain, for e,£ 0 > 0, 

(K{\Vp\)Vp- K{\Vp h \)Vp h ,\e h \ q ~ 2 Vd) <C(\V X \,\9 h \ q - 2 \Vd\) 

< e 0 {K{O |V X | 2 , \9 h \ q - 2 ) + Ce^iK-'&m 2 , \9 h \ q ~ 2 ) 

< eo{K{0\V X \ 2 , \9 h \ q - 2 ) + e \\0 h \\ q Lq + Ce^ ||iWh£)W 

From (|5.15[) , (I5.16P with £o = (q ~ l)(/3 — l)/2, we find that 


Q 

l* 


~ IIMb + {q ^ 1 ] {m) ivxl 2 , \o h \ q ~ 2 ) 


1 d 
q dt 


< Ce 1 - 9 ||i? t ||l, + 2s \\0 h \\l t + Ce~K~^{£)Vd 


(5.15) 


(5.16) 


Q 

L* 


(5.17) 
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Integrating in time from 0 to T 

1 


m 


L°°(I,LQ) 


+ {q 1)( f 1} ^ V(OIVx| 2 , \o h \ q - 2 )dt 


< 2eT Ph\\ q L o. {IM) + Ce 1-9 j ||0 t ||l, dt + Ce~^ J |if-*(£)V0 


dt. 


Li 


Choosing s = using (12.31) we find that 


\\ 0 h\\ q L ^ Lq) + J\m)\v x \ 2 , \e h \ q ~ 2 )ds 


< CT q ~ l / \\d t \\ q Lg dt + CT— / Jr*(f)W 


dt. 


Li 


Note that 


A'-=(£)Vd 


dt = / K~i(C)\X/d\ q dtdx 

Lq Jn J o 



1/2 


Moreover 1 + < C(1 + max{|Vp| a «, \Vp h \ aq })- Thus 

m q L ^ Lq) + f\im\s7 X \\\e h \ q - 2 )dt<CT q ~ l [ T \\Mhdt 

Jo Jo 


q — 2 

+ CT— 


(J J 1 + max{\Vp\ aq ,\Vp h \ aq }dxdty \\Vd\\ 


/o Jn 

Dropping the second term of the LHS in above estimate we obtain 

nl- 


, 1 11 


L 2q(/ ?L 2q) • 


(5.18) 


11^^ IIl c!0 (/,L 9 ) — T T 1 ^9) + CAl T 2 1 || ^d\\ L 2q^J L 2q^ ■ 

Inequality (15.111) follows thanks to (15.181) and triangle inequality 
This finishes the proof. □ 

For the completeness we derive L°°-estimate for p — ph- 

Theorem 5.3. Let p,ph be solutions to problems (E2D and m respectively. Assume 
that p £ L 00 (I 1 W r,00 (Ll)), p t £ L 2 (I,W r ’ q (fl)). Then there exist constant positive constant C 
independence of h such that 


I P Ph\\L°°(I,L°°(Q)) ^ Ch r ||p|lioo(J jly r, 00 (H)) + CT ih q \\Pt\\ L v(p W r, q (Q)} 


+ CA1T 2 ih r ■< \\h'\\L 2 <i(i,w r ’ 29(a)) • 

Proof. Recall that the quasi-uniform of Th we have the inverse estimate (see [HISS]') 

II^IIl“(0) — Ch 1 ||07i|| i( ,(f2) • 

This and triangle inequality imply that 

llxlll/~(/,L°o ) - II^IIl""(/,L-) + 

<Ch ||p|li«j(/ i w r >°°) + Ch q II^^IIl oo (/,l < !) ■ 

Thus inequality (15.191) follows directly from (l5.18l) . □ 


(5.19) 
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5.2. Error estimate for gradient. Now we give an error estimate for gradient 
Theorem 5.4. Under the assumption of Theorem \5.1\ For any 0 < t 0 < t < T. There 
exist constant positive constants C independence of h such that 


\\^{P~ Ph)(t)\\ Lf , {n) < CVM^h ^ \m^ x 


T*h^\\p t {t)\\ 


l 2 (i,h t ( n)) 


} 


+ CM^h L ^~ ||p(f)|| 


(5.20) 


W r <P(Q) 


where the positive constantt M is defined as in (14.171) , and 


V = 


ifo 1 + + J h{t) (l + ||/t(i)|| 2 ) dt 


(5.21) 


Proof. Choosing Wh = 9h in (15.31) we have 


(K(\Vp\)Vp - K(\Vp h \)Vp h ,Vp - S7 Ph ) 

= { Pt - p hit , 0 h ) + (K(\Vp\)Vp - K{\Vp h \)Vp h , VO) . (5.22) 

By (12.91) and (15.51) . there is a positive constant Co independence of h such that 

(K(\Vp\)Vp - K(\Vp h \)Vp h , Vp - Vp h ) > C q M^ ||V(p - p h )f L , . 

In virtue of triangle inequality, Holder’s inequality we have 

(Pt - Ph,u Oh) + (K(\Vp\)Vp - K(\Vp h \)Vp h ,V0) 

< (|/o*| + \Ph,tl\0 h \) + C^Vpf- 1 + \Vp h f~\ |V0|) (5.23) 

< (llPtll + IIpmII) l|0fc|| + CiM?* ||W|| il9 

for some positive constant C\ independence of h. 

Thus 

||V(p - p h )Wl, < C^M^\\p t \\ + ||p M ||) 110,11 + ||W|| L „. 

Thanks to (14.231) and (14.391) . 

INI + \\Ph,t\\ < 22? 2 . (5.24) 


Hence there is a positive constant C independence of h such that 

\\V( P -p h )\\ L , <CVM* \\0 h p + CM^h \\Vd\\l, . (5.25) 

Due to (15.101) and the fact that || Vd|| i/3 < Ch r ~ 1 ||p|| p the left hand side of (I5.25[) is bounded 
by 


CVM 2 


M 2 ' 3X h 2 \\p\\ 2 L 2^ IW r,fi^ + T 2 h r \\pt\\ L 2( IH r) 


+ Ch 2 M 2 ? x 


i w r -f> 


(5.26) 


Thus (15.201) follows from (15.251) and (15.261) . □ 

5.3. Error analysis for fully discrete Galerkin method. In analyzing this method, 
proceed in a similar fashion as for the semidiscrete method, we derive a error estimate for the 
fully discrete time Galerkin approximation the differential equation. 

Let p n {-) = p(-,t n ), be the true solution evaluated at the discrete time levels. We will also 
denote tt p n £ Wh to be the projections of the true solutions at the discrete time levels. As in 
the semidiscrete case, we use x = P~ Ph, 0 = p — Ph, Oh = Ph~^P and \ n , d n , 0% be evaluating 
X, 0, Oh at the discrete time levels. We also define 


d<f> n 


f 1 - cj) n ~ l 


At 
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We rewrite (13.21) with t = t n . Using the definition of L 2 -projection and standard manipulations 
show that the true solution satisfies the discrete equation 

(<9tt p n ,w h ) + (K(\Vp n \)Vp n ,Vw h ) = (- 7 ip™ + dnp n ,w h ) - {if n ,w h ) + {f n ,w h ),\/w h <E W h . 

(5.27) 

Theorem 5.5. Let p n solve problem m and psolve the fully discrete Galerkin finite 
element approximation (ESI) for each time step n, n = 1 There exists a positive 

constant C independent of h and At such that if the At is sufficiently small then 

max ||p" — ph\\ < c (hT^ + V~At) . (5.28) 

l<n<N \ / 


Proof. Subtracting (15.271) from (13.51) gives the error equation 

(dd n h ,w h ) + (if(|Vp£|)Vp£ - tf(|Vp"|)Vp", Vw h ) = (np? - dirp n ,w h ). (5.29) 


Selecting Wh = 0%, we rewrite above equation as form 


+ (K(\Vp n \)Vp n - K(\Vpfl)VpZ,Vp n - S7pl) 

= ( K(\X7p n \)Vp n - K(\Vp n h \)\7ptVd n ) + (tt p? dnp n ,0 n h ). 


For the first term, we have the identity 

( I nn I fj n ~1 Af \ 

do n h , h 2 h + —den) 

= ^ t {n 2 -\K- 1 \\ 2 ) + ^\\d 0 nf- 

For the second term, the monotonicity of K(-) in (12.91) yields 

(lF(|Vp"|)Vp" - K{\S7pl\)Vpl, Vp n - V P £) > C oW " ||V X n ||L 


(5.30) 


(5.31) 


(5.32) 


with uj n = w(t„) = (l + max{||Vp£||^ (n) , II Vp”|| i(5(n) |) . 

For third term, using (12.41) . Holder’s inequality, there is a positive constant C independence of 
h, At, n such that 


(K(\S7p n \)Vp n - K(\V P n\)S7pl,Ve n ) < CM A ||W"|| ifl ■ 
For the last term, it follows from using L 2 -projection and Taylor expand that 

Ptt{ T )( T ~ t n -i)dr, On 

Differentiating equation (11.151) in time provides, 

Ptt = V-{K{\Vp\)Vp) t + f t . 

This and (15.341) imply 

(TTp? - dwp n ,0ft = ^ f t {T)(T - tn-!)dT,0Z 

(. K(\Vp\)Vp) t (r-t n . 




(*p?-dnp n ,ffH) = {p?-dp 



(5.33) 


(5.34) 


(5.35) 
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For the first term of left hand side in (15.3511 , 


h < 


At 




'in-1 


\K 


- At WM\l 2 (i„,l 2 ) 

<iAt||/ t || 2 L2(/nii2) + ^|K|| 2 , 


(r - t n -i) 2 dr 


Wh 


(5.36) 


where/„ = [t„_i,t„]. 

For the second term of left hand side in (15.351) . using integration by part, triangle inequality, 
(EH) , and Holder’s inequality we find that 


12 ^ “7- 

- At 


K(\Vp\)VpdT,VffZ 


tn-1 


\{K{\v P n \)Vp n yoi)\ 


< 


< 


c_ 

At 

C 
A t 


tn-1 


(\\7p\ 0 ~\ |V0£|) dr + C(\Vp n f~\ |V^|) 


>t n - 1 


IIVpHi^dr 


l|V^||^+C'||Vp™||^ 1 ||V^|| L ,. 


According to (14.381) . 


/ 2 < GM~p~ ||V0£|| iP 


< CM 2 ^ || Vx”II lp +CM V ||W"|Il/» 


(5.37) 


< e || V X n || 2 ^ + Ce-'M 2 ^ 1 + CM^ ||W"|| L * . 

Combining (15.351) . (15.361) and (15.371) and using the fact that j = 

(np?-dnp n ,eZ) < ^Ai ||/t||i 2 ( j ni x, 3) + \ m 2 

+ £ ||Vx"||L> + Ce^Mi + CM * ||V0"|Il/> . 

It follows from (15.311) . (15.321) . (15.331) and (15.381) that 

2 ^ (ll^ll 2 - ll^" 1 !! 2 ) + CW" IIVxll* < £ l|Vx"||^ + C(M A + Mi) ||V^|| L , 

+CAt ||/t||i2(j n|i 2) + 2 Ill'll 2 + Ge l M^. 

Selecting e = c '°^ in (15.391) . we obtain 

^ (ll^ll 2 - II^M 2 ) + 0)W n l|Vx n ||^ < C(M^ + Mi) ||Vt9"|| L/) 40) 

+ CAt + ||^|| 2 + Ciuj^M^ + CMi || V r ?"|| L p . 

Dropping the second term in (15.401) and using (15.51) and the fact j- + 7 = 1, we get 

^ (IKH 2 - IK _1 H 2 ) < ll^ll 2 + C {{M^ + Mi) ||V^|L, + X l+7 } + CAt ||/i|| 2 i2(7 „, L2) ■ 
Summing from n = 0 to n = m—1 with 0 < m < N and using the fact that 0^(0) = 0, we 


(5.38) 


(5.39) 


obtain 


n—0 


Wf < A \K\\ 2 + c{(M^+Mi) At ||W "|| L 0 

n —0 

m— 1 

+ (At ) 2 ll/tl| 2 2 (/„,L 2 ) + {m- l)A4At|. 


n —0 


(5.41) 
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An application of the discrete Gronwall’s inequality shows that for At sufficiently small, 


m —1 


o maxjra 2 <c{(M^+M$) £ At ||W"||^ + MAt} + C(Atf \\f t \\ 2 L2(I L2) , 


n —0 


which implies 


(5.42) 


iV-1 1 

max^ ||0"|| < c{(M w + M™)br? ( £ At \\p n \\ wr ,, ) 3 +M*VAt} 

- n ~ V n—0 (5.4«5j 

+ ||/t|| i 2(J ii 2) • 

Finally, the result follows from (15.431) and \\p n — pJJ|| < W®h\\ + ||d"|| . □ 

Now we give the error estimate for gradient in fully discrete version. 

Theorem 5.6. Let p solve problem m and p^ solve the fully discrete Galerkin finite 
element approximation (13751) for each time step n, n = 1,..., N. Then there exists a positive 
constant C independent of h and At satisfying 

II - Vp n \\ LP < C(h ^ + VAt). (5.44) 

Proof. Selecting Wh = 6%, in (15.291) and rewriting the resulting equation as form 

((tf(| Vp"|)Vp" - K(\Vpft)Vph V/U - S7p n h ) 

= (K(\S7 p n \)S7 p n - K{\Vp n h \)Vp n h , W") + (p? - dp n h , 6%). 

Due to (12.101) . Cauchy-Schwartz and triangle inequality, one has 

w" II- vp"II 2 , < c (\\dpi\\ + M\\) we n h \\ + c{ |Vp" - |W"|). 


Since 


ll» = (At) -1 


Ph,tdt 


< (A t )- 1 


\\Ph,t\\ dt < max \\ph,t\\ < P 2 , 

[T/N,T] " 


||p"11 < max[ T /jv,T] Ibt|| < P 2 and (15.331) . we obtain 

oc" II Vpl - Vp n Wl, < CP 2 ||6>J*|| +CM A HWli, 
Using estimate in (15.431) and (15.51) shows that 


JV-l i 

\\VPh ~ Vp"!! 2 , < CP 2 M^{(M^ + M^)h^( At \\p n \\wr,e ) 5 + m^VaI 

71=0 

+ At ||/t|| L 2 (Jii3) } + CM^hT - 1 Wp n \\w^ : 

which proves (15.441) . The proof is complete. □ 

6. Numerical results. In this section, we give a two simple numerical experiments using 
Galerkin finite element method in the two dimensional region to illustrate the convergent 
theory. For simplicity, we test the convergence of our method with the Forchheimer two-term 
law g(s) = 1 + s. Equation (11.121) sg(s) = £, s > 0 gives s = ~ 1+ ^ 1+4 - and hence 


g(s(f)) 1 + v /T+4T 

The region is selected is unit square, i.e U = [0, l] 2 . We use the Lagrange element of order 
r = 2 on the unit square in two dimensions. We used FEniCS 1277 to perform our numerical 
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simulations. We divide the unit square into an N x N mesh of squares, each then subdivide 
into two right triangles using the UnitSquareMesh class in FEniCS. For each mesh, we solve 
the generalized Forchheimer equation numerically. The error control in each nonlinear solve 
is e = 10~ 6 . Our problem is solved at each time level start at t = 0 until final time T = 1. 
At this time, we measured the L 2 -errors of pressure and L^-errors of gradient of pressure and 
velocity. Here = 2- a = 2- d tg(g)h = §• 

Example 1. The analytical solution is as follows 


p(x,t) = e 24 





V(x, i)eiix [o, l]. 


The forcing term / is determined accordingly to the analytical solution by equation p t — V • 
(K(\S7p\)p) = f. Explicitly, 


f(x,t) = —2e 24 ^(xl +x%) - i(x? +x\) 


4e 24 (1 — x\ — X 2 ) 


1+ v/1 + 4 e~ 2t W(x) 

4e~ 4t \x\ (1 — xi) 2 (l — 2xi) + x 2 (l — x 2 ) 2 (l — 2x 2 )l 


W{x) (l + yjl + 4e~ 2t W(xf) y/1 + 4e- 2t W(x) 


where W(x) = \Zxf(l — xi) 2 + x|(l — x 2 ) 2 . The initial data p°(x) = |(x 2 +X 2 ) — ^(xf — x^ + l 
and the Neumann boundary condition ip(x,t) = 0. The numerical results are listed as the 
following table. 


N 

IIP-Phil 

Rates 

ll^(P Ph)|| lP(Q) 

Rates 

4 

6.33e-02 

- 

4.51E-01 

- 

8 

5.50E-02 

0.20 

4.07E-01 

0.15 

16 

4.52E-02 

0.28 

3.70E-01 

0.14 

32 

3.50E-02 

0.37 

3.22E-01 

0.20 

64 

2.53E-02 

0.47 

2.70E-01 

0.25 

128 

1.73E-02 

0.55 

2.21E-01 

0.29 

256 

1.13E-02 

0.61 

1.79E-01 

0.31 

512 

7.19E-03 

0.63 

1.43E-01 

0.32 


Table 1. Convergence study for generalized Forchheimer flows using Galerkin finite element 
method with zero flux on the boundary in 2D. 

Example 2. The analytical solution is p(x, t) = xix 2 e -4 + 1 for all (x, t) G D x [0, 1]. The 
forcing term /, initial condition and Neumann boundary condition are determined accordingly 
to the analytical solution as follows 


/(x, t) = —e t X\X 2 + 


8 e 2t X\X 2 


\Jx\+ x\ ^1 + \j\ + 4e l \/x\ + xfj \J 1 + 4e 4 \/x\ + 


x% 


and 


p°(x)=x ix 2 + l, ip(x,t) = 


2 e~ 


1 + \J 1 + 4 e~ t ^Jx\ + ; 


x 2 

on xi 

= 0 

-X2 

on x\ 

= 1 

-Xi 

on #2 

= 1 

Xi 

on #2 

= 0 


The numerical results are listed in below table 
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N 

\\P~ Ph\\ 

Rates 


Rates 

4 

4.40E-02 

- 

2.67E-02 

- 

8 

2.24E-02 

0.97 

2.02E-02 

0.40 

16 

1.15E-02 

0.96 

1.37E-02 

0.56 

32 

5.90E-03 

0.96 

8.53E-03 

0.68 

64 

3.01E-03 

0.97 

4.99E-03 

0.77 

128 

1.53E-03 

0.98 

2.79E-03 

0.84 

256 

7.70E-04 

0.99 

1.52E-03 

0.88 

512 

3.94E-4 

0.99 

8.28E-4 

0.92 


Table 2. Convergence study for generalized Forchheimer flow using Galerkin finite element 
method with nonzero flux on the boundary in 2D. 

Acknowledgment. The author wishes to thank Dr. Luan Hoang for valuable assistance 
in discussion and suggestions. 


REFERENCES 

[1] T. Arbogast, M. F. Wheeler, and N.-Y. Zhang, A nonlinear mixed finite element method for a 

degenerate parabolic equation arising in flow in porous media , SIAM J. Numer. Anal., 33 (1996), 
pp. 1669-1687. 

[2] E. Aulisa, L. Bloshanskaya, L. Hoang, and A. Ibragimov, Analysis of generalized Forchheimer flows 

of compressible fluids in porous media , J. Math. Phys., 50 (2009), pp. 103102, 44. 

[3] J. Bear, Dynamics of Fluids in Porous Media , Dover, New York, 1972. 

[4] S. C. Brenner and L. R. Scott, The mathematical theory of finite element methods , vol. 15 of Texts in 

Applied Mathematics, Springer, New York, third ed., 2008. 

[5] F. Brezzi and M. Fortin, Mixed and hybrid finite element methods , vol. 15 of Springer Series in 

Computational Mathematics, Springer-Verlag, New York, 1991. 

[6] P. G. Ciarlet, The finite element method for elliptic problems, North-Holland Publishing Co., 

Amsterdam, 1978. Studies in Mathematics and its Applications, Vol. 4. 

[7] E. DiBenedetto, Partial differential equations, Cornerstones, Birkhauser Boston Inc., Boston, MA, 

second ed., 2010. 

[8] J. J. DOUGLAS, P. J. Paes-Leme, and T. Giorgi, Generalized Forchheimer flow in porous media, in 

Boundary value problems for partial differential equations and applications, vol. 29 of RMA Res. 
Notes Appl. Math., Masson, Paris, 1993, pp. 99-111. 

[9] K. B. Fadimba, Error analysis for a Galerkin finite element method applied to a coupled nonlinear 

degenerate system of advection-diffusion equations, Comput. Methods Appl. Math., 6 (2006), pp. 3-30 
(electronic). 

[10] -, On existence and uniqueness for a coupled system modeling immiscible flow through a porous 

medium, J. Math. Anal. Appl., 328 (2007), pp. 1034-1056. 

[11] -, A linear backward Euler scheme for the saturation equation: regularity results and consistency, 

J. Comput. Appl. Math., 234 (2010), pp. 272-282. 

[12] -, A linear backward Euler scheme for a class of degenerate advection-diffusion equations: a 

mathematical analysis of the convergence in L°°(0, To; L 2 (f2)) and in L 2 (0, To; H 1 (£1)), Anal. Appl. 
(Singap.), 12 (2014), pp. 227-249. 

[13] K. B. Fadimba and R. C. Sharpley, A priori estimates and regularization for a class of porous medium 

equations, Nonlinear World, 2 (1995), pp. 13-41. 

[14] -, Galerkin finite element method for a class of porous medium equations, Nonlinear Anal. Real 

World Appl., 5 (2004), pp. 355-387. 

[15] P. Forchheimer, Wasserbewegung, Zeit. Ver. Deut. Ing., 45 (1901), pp. 1781-1788. 

[16] P. Forchheimer, Hydraulik, no. Leipzig, Berlin, B. G. Teubner, 1930. 3rd edition. 

[17] L. Hoang and A. Ibragimov, Structural stability of generalized Forchheimer equations for compressible 

fluids in porous media, Nonlinearity, 24 (2011), pp. 1-41. 

[18] L. Hoang and A. Ibragimov, Qualitative Study of Generalized Forchheimer Flows with the Flux 

Boundary Condition, Adv. Diff. Eq., 17 (2012), pp. 511-556. 

[19] L. Hoang, A. Ibragimov, T. Kieu, and Z. Sobol, Stability of solutions to generalized Forchheimer 

equations of any degree, J. Math. Sci., (2015). accepted. 

[20] L. Hoang and T. Kieu, Global estimates for generalized Forchheimer flows of slightly compressible fluids, 

Journal d’Analyse Mathematique, (2015). accepted. 

[21] -, Interior estimates for generalized Forchheimer flows of slightly compressible fluids, (2015). 

submitted, preprint http://arxiv.org/abs/1404.6517 

[22] L. T. Hoang, A. Ibragimov, and T. T. Kieu, One-dimensional two-phase generalized Forchheimer flows 

of incompressible fluids, J. Math. Anal. Appl., 401 (2013), pp. 921-938. 

[23] -, A family of steady two-phase generalized Forchheimer flows and their linear stability analysis, J. 

Math. Phys., 55 (2014), p. 123101. 















Galerkin finite element method for Forchheimer equation of slightly compressible fluids 


23 


[24] L. T. Hoang, T. T. Kieu, and T. V. Phan, Properties of generalized Forchheimer flows in porous media, 

J. Math. Sci., 202 (2014), pp. 259-332. 

[25] C. Johnson and V. Thomee, Error estimates for some mixed finite element methods for parabolic type 

problems , RAIRO Anal. Numer., 15 (1981), pp. 41-78. 

[26] T. Kieu, Analysis of expanded mixed finite element methods for the generalized Forchheimer flows of 

slightly compressible fluids, Numer. Methods Partial Differential Equations, (2015). accepted. 

[27] O. A. Ladyzenskaja, V. A. Solonnikov, and N. N. Uralceva, Linear and quasilinear equations of 

parabolic type , Translated from the Russian by S. Smith. Translations of Mathematical Monographs, 
Vol. 23, American Mathematical Society, Providence, R.I., 1968. 

[28] J.-L. Lions, Quelques methodes de resolution des problemes aux limites non lineaires, Dunod, 1969. 

[29] A. Logg, K.-A. Mardal, and G. N. Wells, eds., Automated Solution of Differential Equations by the 

Finite Element Method, vol. 84 of Lecture Notes in Computational Science and Engineering, Springer, 

2012 . 

[30] M. Muskat, The flow of homogeneous fluids through porous media, McGraw-Hill Book Company, inc., 

1937. 

[31] D. A. Nield and A. Bejan, Convection in porous media, Springer-Verlag, New York, second ed., 1999. 

[32] E.-J. Park, Mixed finite element methods for generalized Forchheimer flow in porous media, Numer. 

Methods Partial Differential Equations, 21 (2005), pp. 213-228. 

[33] R. E. Showalter, Monotone operators in Banach space and nonlinear partial differential equations, 

vol. 49 of Mathematical Surveys and Monographs, American Mathematical Society, Providence, RI, 
1997. 

[34] B. Straughan, Stability and wave motion in porous media, vol. 165 of Applied Mathematical Sciences, 

Springer, New York, 2008. 

[35] V. Thomee, Galerkin finite element methods for parabolic problems, vol. 25 of Springer Series in 

Computational Mathematics, Springer-Verlag, Berlin, second ed., 2006. 

[36] J. C. Ward, Turbulent flow in porous media., Journal of the Hydraulics Division, Proc. Am. Soc. Civ. 

Eng., 90(HY5) (1964), pp. 1-12. 

[37] C. S. Woodward and C. N. Dawson, Analysis of expanded mixed finite element methods for a nonlinear 

parabolic equation modeling flow into variably saturated porous media, SIAM J. Numer. Anal., 37 
(2000), pp. 701-724 (electronic). 

[38] E. Zeidler, Nonlinear functional analysis and its applications. II/B, Springer-Verlag, New York, 1990. 

Nonlinear monotone operators, Translated from the German by the author and Leo F. Boron. 



